/**
我们将会使用一个例子来进行CPU与GPU之间的异步计算，并且进行适当地同步工作。我们将利用GPU计算两个数组的大小，
分别使用流来进行并行计算。在每一个block中我们将利用共享内存进行规约求和，不同block之间我们将使用全局内存进行并行串行求和。
在这个过程中，不仅会用到共享内存来进行规约求和，也会用到共享变量来记录当前已经完成的求和情况，方便最后对每一个block结果的串行归约。
我们需要根据输出数组的情况开辟合适的线程网格来解决问题。并且使用事件来进行时间的记录。这个例子将会利用到几乎所有的CUDA的基本内容。
**/


#include<stdio.h>
#include<stdlib.h>
#include<string.h>

//需要求和的数组大小
#define N 1000

//用来记录计算完成情况的变量
__device__ int count1 = 0;
__device__ int count2 = 0;

//判断是不是可以将过程量相加的变量
// __device__ int judge = 0;

//首先声明内核函数，这个函数处理的就是进行一个数组中所有元素的相加操作
//传入参数就是当前数组的指针以及当前数组的大小。
//并且传入一个数组result来保存每个block计算的中间结果，将最后的计算结果存在这个中间结果数组的第0号位置。
//中间结果数组的大小为gridDim.x
//第四个参数传入流的编号
__global__ void arr_sum(int* input_arr, int size, long* result, int stream_id){
    // printf("进入内核函数\n");
    //这里进入内核函数
    //首先先进行规约求和，首先先申请共享内存，一共192位
    //注意关键字是shared而不是share
    __shared__ long temp_arr[192];

    //计算总的线程编号，来索引全局内存
    unsigned int thread_idx = threadIdx.x + blockIdx.x * blockDim.x;

    // //进行数据拷贝，因为数组的大小可能不是对齐的，加了判断防止数组越界
    if(thread_idx < size){
        temp_arr[threadIdx.x] = input_arr[thread_idx];
    }else {
        //因为输出的数组可能
        temp_arr[threadIdx.x] = 0;
    }

    // printf("共享内存拷贝完毕\n");

    // 拷贝结束之后必须进行同步，保证所有warp的操作全部完成
    __syncthreads();

    //这里进行规约求和，规约求和的特点就是每个线程只执行一个计算
    //所以只有数据规模一半的线程是活跃的
    int active = 192 / 2;

    //这里开始进行规约求和
    //我们需要一个循环，让工作线程的边界不断缩减
    int i;
    //当规约到奇数个项的时候就难以规约了，这个时候使用暴力相加的方式就好
    //当active的边界是3的时候，还有6个元素，还可以规约
    //所有i的边界包括3
    for(i = active; i >= 3; i = i / 2){
        //如果当前线程是活跃线程就参与计算
        if(threadIdx.x < i){
            //开始计算两个点的相加
            temp_arr[threadIdx.x] = temp_arr[threadIdx.x] + temp_arr[threadIdx.x + i];
        }
        
        //不要在分支中同步
        __syncthreads();
    }

    // // printf("规约求和完毕\n");
    
    int judge = 0;
    if(threadIdx.x == 0){
        //然后使用每个block的第一个线程来做最后三个元素的相加
        // printf("%ld,%ld,%ld\n", temp_arr[0], temp_arr[1], temp_arr[2]);
        result[blockIdx.x] = temp_arr[0] + temp_arr[1] + temp_arr[2];
        // printf("%ld\n", result[blockIdx.x]);
        //这里强制等待写回
        __threadfence();
        //这里记录一个变量，每个块完成计算之后，就为这个值自增
        if(stream_id == 1){
            atomicAdd(&count1,1);
            if(count1 == gridDim.x){
                judge = 1;
            }
        }else{
            atomicAdd(&count2,1);
            if(count2 == gridDim.x){
                judge = 1;
            }
        }
    }

    //这里代表最后一个块，最后一个块进行收尾
    if(judge == 1){
        if(threadIdx.x == 0){
            for(i = 1; i < gridDim.x; i++){
                // printf("%ld\n", result[i]);
                result[0] = result[0] + result[i];
            }
        }
    }
    
}

//主函数
int main(){
    //这里计算网格结构，我们使用的是Tesla K40，15SM，192SP
    //我们建立一维网格，计算需要几个block
    int block_num = (N / 192);

    //如果数组的大小没有办法与192取整，那么就需要一个新的块
    if(N % 192 != 0){
        block_num++;
    }
    
    //创建两个数组为了方便验证结果，分别是正序的和倒序的等差数列。数组的大小为N
    //使用pinned memory
    int* arr1;
    int* arr2;
    long* result1;
    long* result2;
    cudaMallocHost((void **)&arr1, N * sizeof(int), cudaHostAllocWriteCombined);
    cudaMallocHost((void **)&arr2, N * sizeof(int), cudaHostAllocWriteCombined);
    //用来存储结果的数组
    cudaMallocHost((void **)&result1, block_num * sizeof(long));
    cudaMallocHost((void **)&result2, block_num * sizeof(long));
    

    //用来迭代的变量
    int i;
    //将两个数组分别初始化
    for(i = 0; i < N; i++){
        arr1[i] = i;
        arr2[i] = N - i - 1;
    }

    printf("主机数据初始化完毕\n");


    //申请GPU空间
    int* device_input_arr1 = NULL;
    int* device_input_arr2 = NULL;
    long* device_result1 = NULL;
    long* device_result2 = NULL;


    //使用cudaMalloc来申请空间
    //函数原型__host____device__cudaError_t cudaMalloc (void** devPtr, size_t size)
    //注意这里要传入的是指针的地址
    cudaMalloc((void **)&device_input_arr1, N*sizeof(int));
    cudaMalloc((void **)&device_input_arr2, N*sizeof(int));
    cudaMalloc((void **)&device_result1, block_num*sizeof(long));
    cudaMalloc((void **)&device_result2, block_num*sizeof(long));

    printf("设备空间分配完毕\n");

    //创建两个流
    cudaStream_t stream[2];

    for(i = 0; i < 2; i++){
        cudaStreamCreate(&stream[i]);
    }

    //完成数据拷贝，从这里开始使用流
    cudaMemcpyAsync((void *)device_input_arr1, (void *)arr1, N*sizeof(int), cudaMemcpyHostToDevice, stream[0]);
    cudaMemcpyAsync((void *)device_input_arr2, (void *)arr2, N*sizeof(int), cudaMemcpyHostToDevice, stream[1]);
    
    //初始化中间结果数组
    cudaMemsetAsync((void *)device_result1, 0, block_num * sizeof(long), stream[0]);
    cudaMemsetAsync((void *)device_result2, 0, block_num * sizeof(long), stream[1]);

    printf("激活初始化\n");
    //网格建立，一维网格
    dim3 grid_dimension(block_num);
    //每个block拥有和sp数量一致的thread，也是一维的
    dim3 block_dimension(192);

    cudaError_t  error_check;

    //使用内核函数，不分配动态共享内存
    arr_sum<<<grid_dimension, block_dimension, 0, stream[0]>>>(device_input_arr1, N, device_result1, 1);
    arr_sum<<<grid_dimension, block_dimension, 0, stream[1]>>>(device_input_arr2, N, device_result2, 2);

    error_check = cudaGetLastError();
    if( error_check != cudaSuccess ){
        printf("%s\n" , cudaGetErrorString( error_check ) );
        system("pause") ;
        return 0 ;
    }

    printf("激活内核\n");
    //将数据拷贝出来
    cudaMemcpyAsync((void *)result1, (void *)device_result1, block_num * sizeof(long), cudaMemcpyDeviceToHost, stream[0]);
    cudaMemcpyAsync((void *)result2, (void *)device_result2, block_num * sizeof(long), cudaMemcpyDeviceToHost, stream[1]);
    printf("激活数据拷出\n");
    printf("进入等待\n");
    //这里等待所有的GPU代码执行完毕
    cudaThreadSynchronize();

    //检查数据是不是正确
    long real_result = (0 + N - 1) * N / 2;

    printf("结果%ld\n", result2[0]);

    if(result1[0] == real_result){
        printf("第一个结果正确\n");
    }else{
        printf("第一个结果不正确\n");
    }

    if(result2[0] == real_result){
        printf("第二个结果正确\n");
    }else{
        printf("第二个结果不正确\n");
    }
    

    //析构两个流
    for(i = 0; i < 2; i++){
        //注意这个函数的使用，流的创建输入的是
        cudaStreamDestroy(stream[i]);
    }

    //析构pinned memory
    cudaFreeHost(arr1);
    cudaFreeHost(arr2);
    cudaFreeHost(result1);
    cudaFreeHost(result2);

    //析构在设备端分配的函数
    cudaFree(device_input_arr1);
    cudaFree(device_input_arr2);
    cudaFree(device_result1);
    cudaFree(device_result2);

    return 0;
}